;;..............................................................................  
;; merge U,V,T,Z3 and calculate and S     
;;..............................................................................  
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

   yy = "2010" 
   mm = "04" 
   res = "0.9x1.25" 
   res = "1.9x2.5" 
   res = "2x2.5" 

   load "loadit" 

   tpath = "/glade/scratch/kaizhang/data/ndata/tmp/" 
   ipath = "/glade/scratch/kaizhang/data/ndata/input/" 
   opath = "/glade/scratch/kaizhang/data/ndata/out/" 
   spath = "/glade/scratch/kaizhang/data/ndata/cam5nud/" 

   fna = opath+"ndg_v1_"+res+"_"+yy+mm+".nc" 
   fnb = opath+"z3_"+res+"_" +yy+mm+".nc" 
   fnc = opath+"ecmwf_1.9x2.5L30_200901_ori.nc" 
   fno = opath+"ecmwf_1.9x2.5L30_"+yy+mm+".nc" 
  
   print("input files : "+fna) 
   ;;print("input files : "+fnb) 
   print("input files : "+fnc) 
 
   fla = addfile(fna,"r") 
   flb = addfile(fnb,"r") 
   flc = addfile(fnc,"r") 

   flo = addfile(fno,"w") 

;;..............................................................................  
;; read data 
;;..............................................................................  
   print("get ECMWF data : T") 
   vt = fla->var130 ;;T 
   print("get ECMWF data : U") 
   vu = fla->var131 ;;U 
   print("get ECMWF data : V") 
   vv = fla->var132 ;;V
   print("get ECMWF data : PS") 
   vp = fla->var152 ;;PS
   hyam = flc->hyam
   hybm = flc->hybm
   P0   = flc->P0
   vzs  = fla->var156 ;;flb->Z3
   vz   = flb->Z3

   nt = dimsizes(vt(:,0,0,0)) 

;;..............................................................................  
;; Get cam5 nudging data frame 
;;..............................................................................  
   print("get nudging frame data : T") 
   T = flo->T 
   print("get nudging frame data : U") 
   U = flo->U 
   print("get nudging frame data : V") 
   V = flo->V
   print("get nudging frame data : S") 
   S = flo->S
   print("get nudging frame data : PS") 
   PS= flo->PS 

;;..............................................................................  
;; ECMWF data 
;;..............................................................................  
   T(0:nt-1,:,:,:) = (/ flt2dble(vt(:,:,:,:)) /) 
   U(0:nt-1,:,:,:) = (/ flt2dble(vu(:,:,:,:)) /) 
   V(0:nt-1,:,:,:) = (/ flt2dble(vv(:,:,:,:)) /) 
  PS(0:nt-1,:,:) = (/ flt2dble(vp(:,:,:)) /) 

   g  = 9.8d0     ;; m s-2
   cp = 1005.7d0  ;; J kg-1 K-1

   S = T
   S@long_name = "Dry static energe"
   S@units     = "m2 s-2"
   ;;S = (/ cp*T + g*Z3 /)

   printVarSummary(vt) 
   printVarSummary(vz) 
   printVarSummary(S) 

   nlev = 30 

;;..............................................................................  
;; Compute S, remember to add surface geopotential 
;;..............................................................................  
   do ik = 0,nlev-1 
      S(0:nt-1,ik,:,:) = (/ cp * flt2dble(vt(:,ik,:,:)) + g*flt2dble(vz(:,ik,:,:)+vzs(:,:,:)) /) 
   end do 

;;..............................................................................  
;; additional information 
;;..............................................................................  
   T@source = "ECMWF operational analysis" 
   U@source = "ECMWF operational analysis" 
   V@source = "ECMWF operational analysis" 
   S@source = "ECMWF operational analysis" 
   PS@source = "ECMWF operational analysis" 

;;..............................................................................  
;; Output data 
;;..............................................................................  
   print("output the data ") 
   flo->T=T 
   flo->U=U 
   flo->V=V 
   flo->S=S 
   flo->PS=PS 
   flo->hyam=hyam
   flo->hybm=hybm
   flo->P0=P0

end





